function U = poisson_solver_2d(mesh, pde)
A = assemble_matrix_2d(mesh, pde);
b = assemble_vector_2d(mesh, pde);
[A, b] = process_boundary_2d(A, b, mesh, pde);
U = A\b;
end

%% assemble_matrix_2d
function matrix_A = assemble_matrix_2d(mesh, pde)
matrix_A = sparse(mesh.Nb,mesh.Nb);
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(V, "2d");
    phix = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type, "dx");
    phiy = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type, "dy");
    fc = pde.cor(gauss_points_loc(:,1), gauss_points_loc(:,2));
    for i = 1:mesh.Nlb
        for j = 1:mesh.Nlb
            f = fc .* (phix{i} .* phix{j} + phiy{i} .* phiy{j});
            r = gauss_weights_loc * f;
            matrix_A(mesh.elem(n,i), mesh.elem(n,j)) = matrix_A(mesh.elem(n,i), mesh.elem(n,j)) + r;
        end
    end
end
end

%% assemble_vector_2d
function vector_b = assemble_vector_2d(mesh, pde)
vector_b = zeros(mesh.Nb,1);
for n = 1:mesh.Ne
    E = mesh.node(mesh.elem(n,:),:);
    V = E(1:3,:);
    [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(V, "2d");
    phi = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type);
    ff = pde.f(gauss_points_loc(:,1), gauss_points_loc(:,2));
    for j = 1:mesh.Nlb
        f =  ff .* phi{j};
        r = gauss_weights_loc * f;
        vector_b(mesh.elem(n,j),1) = vector_b(mesh.elem(n,j),1) + r;
    end
end
end